######################################################
######################################################
#
#		Methods Paper -- Plots for Outcputs
#
#		Lucas Leemann, ltl2108@columbia.edu
#
######################################################
######################################################



# Make Big Data -- SAMPLE SIZE 1000

rm(list=ls())
setwd("/Users/lleemann/Desktop/Replication Files/Simulation/Resulting data")

###### XXXX_-XX_005

## -0.9
load("ResultSimulation_1000_-09_005.RData")

# b
param.struc_m9.005 <- stuff$param.struc
param.selec_m9.005 <- stuff$param.selec
param.selStruc_m9.005 <- stuff$param.selStruc

# c
correlation.005 <- matrix(NA,1000, 5)
correlation.005[,1] <- stuff$corr.rho


## -0.5
load("ResultSimulation_1000_-05_005.RData")

# b
param.struc_m5.005 <- stuff$param.struc
param.selec_m5.005 <- stuff$param.selec
param.selStruc_m5.005 <- stuff$param.selStruc

# c
correlation.005[,2] <- stuff$corr.rho


## +0.0
load("ResultSimulation_1000_+00_005.RData")

# b
param.struc_p0.005 <- stuff$param.struc
param.selec_p0.005 <- stuff$param.selec
param.selStruc_p0.005 <- stuff$param.selStruc

# c
correlation.005[,3] <- stuff$corr.rho


## +0.5
load("ResultSimulation_1000_+05_005.RData")

# b
param.struc_p5.005 <- stuff$param.struc
param.selec_p5.005 <- stuff$param.selec
param.selStruc_p5.005 <- stuff$param.selStruc

# c
correlation.005[,4] <- stuff$corr.rho


## +0.9
load("ResultSimulation_1000_+09_005.RData")

# b
param.struc_p9.005 <- stuff$param.struc
param.selec_p9.005 <- stuff$param.selec
param.selStruc_p9.005 <- stuff$param.selStruc

# c
correlation.005[,5] <- stuff$corr.rho




###### XXXX_-XX_010

## -0.9
load("ResultSimulation_1000_-09_010.RData")

# b
param.struc_m9.010 <- stuff$param.struc
param.selec_m9.010 <- stuff$param.selec
param.selStruc_m9.010 <- stuff$param.selStruc

# c
correlation.010 <- matrix(NA,1000, 5)
correlation.010[,1] <- stuff$corr.rho


## -0.5
load("ResultSimulation_1000_-05_010.RData")

# b
param.struc_m5.010 <- stuff$param.struc
param.selec_m5.010 <- stuff$param.selec
param.selStruc_m5.010 <- stuff$param.selStruc

# c
correlation.010[,2] <- stuff$corr.rho


## +0.0
load("ResultSimulation_1000_+00_010.RData")

# b
param.struc_p0.010 <- stuff$param.struc
param.selec_p0.010 <- stuff$param.selec
param.selStruc_p0.010 <- stuff$param.selStruc

# c
correlation.010[,3] <- stuff$corr.rho


## +0.5
load("ResultSimulation_1000_+05_010.RData")

# b
param.struc_p5.010 <- stuff$param.struc
param.selec_p5.010 <- stuff$param.selec
param.selStruc_p5.010 <- stuff$param.selStruc

# c
correlation.010[,4] <- stuff$corr.rho


## +0.9
load("ResultSimulation_1000_+09_010.RData")

# b
param.struc_p9.010 <- stuff$param.struc
param.selec_p9.010 <- stuff$param.selec
param.selStruc_p9.010 <- stuff$param.selStruc

# c
correlation.010[,5] <- stuff$corr.rho



###### XXXX_-XX_020

## -0.9
load("ResultSimulation_1000_-09_020.RData")

# b
param.struc_m9.020 <- stuff$param.struc
param.selec_m9.020 <- stuff$param.selec
param.selStruc_m9.020 <- stuff$param.selStruc

# c
correlation.020 <- matrix(NA,1000, 5)
correlation.020[,1] <- stuff$corr.rho


## -0.5
load("ResultSimulation_1000_-05_020.RData")

# b
param.struc_m5.020 <- stuff$param.struc
param.selec_m5.020 <- stuff$param.selec
param.selStruc_m5.020 <- stuff$param.selStruc

# c
correlation.020[,2] <- stuff$corr.rho


## +0.0
load("ResultSimulation_1000_+00_020.RData")

# b
param.struc_p0.020 <- stuff$param.struc
param.selec_p0.020 <- stuff$param.selec
param.selStruc_p0.020 <- stuff$param.selStruc

# c
correlation.020[,3] <- stuff$corr.rho


## +0.5
load("ResultSimulation_1000_+05_020.RData")

# b
param.struc_p5.020 <- stuff$param.struc
param.selec_p5.020 <- stuff$param.selec
param.selStruc_p5.020 <- stuff$param.selStruc

# c
correlation.020[,4] <- stuff$corr.rho


## +0.9
load("ResultSimulation_1000_+09_020.RData")

# b
param.struc_p9.020 <- stuff$param.struc
param.selec_p9.020 <- stuff$param.selec
param.selStruc_p9.020 <- stuff$param.selStruc

# c
correlation.020[,5] <- stuff$corr.rho




###### XXXX_-XX_100

## -0.9
load("ResultSimulation_1000_-09_100.RData")

# b
param.struc_m9.100 <- stuff$param.struc
param.selec_m9.100 <- stuff$param.selec
param.selStruc_m9.100 <- stuff$param.selStruc

# c
correlation.100 <- matrix(NA,1000, 5)
correlation.100[,1] <- stuff$corr.rho


## -0.5
load("ResultSimulation_1000_-05_100.RData")

# b
param.struc_m5.100 <- stuff$param.struc
param.selec_m5.100 <- stuff$param.selec
param.selStruc_m5.100 <- stuff$param.selStruc

# c
correlation.100[,2] <- stuff$corr.rho


## +0.0
load("ResultSimulation_1000_+00_100.RData")

# b
param.struc_p0.100 <- stuff$param.struc
param.selec_p0.100 <- stuff$param.selec
param.selStruc_p0.100 <- stuff$param.selStruc

# c
correlation.100[,3] <- stuff$corr.rho


## +0.5
load("ResultSimulation_1000_+05_100.RData")

# b
param.struc_p5.100 <- stuff$param.struc
param.selec_p5.100 <- stuff$param.selec
param.selStruc_p5.100 <- stuff$param.selStruc

# c
correlation.100[,4] <- stuff$corr.rho


## +0.9
load("ResultSimulation_1000_+09_100.RData")

# b
param.struc_p9.100 <- stuff$param.struc
param.selec_p9.100 <- stuff$param.selec
param.selStruc_p9.100 <- stuff$param.selStruc

# c
correlation.100[,5] <- stuff$corr.rho


#	#	#	#	#	#	#	#	#	#	#	#	#	#	#	#	
#	#	Build large things:
#	#	#	#	#	#	#	#	#	#	#	#	#	#	#	#	


#		STRUCTURAL
Param.struc.1000 <- array(NA,c(1000,6,20))
# -0.9 corr
Param.struc.1000[,,1] <- param.struc_m9.005
Param.struc.1000[,,2] <- param.struc_m9.010
Param.struc.1000[,,3] <- param.struc_m9.020
Param.struc.1000[,,4] <- param.struc_m9.100
# -0.5 corr
Param.struc.1000[,,5] <- param.struc_m5.005
Param.struc.1000[,,6] <- param.struc_m5.010
Param.struc.1000[,,7] <- param.struc_m5.020
Param.struc.1000[,,8] <- param.struc_m5.100
# +0.0 corr
Param.struc.1000[,,9] <- param.struc_p0.005
Param.struc.1000[,,10] <- param.struc_p0.010
Param.struc.1000[,,11] <- param.struc_p0.020
Param.struc.1000[,,12] <- param.struc_p0.100
# +0.5 corr
Param.struc.1000[,,13] <- param.struc_p5.005
Param.struc.1000[,,14] <- param.struc_p5.010
Param.struc.1000[,,15] <- param.struc_p5.020
Param.struc.1000[,,16] <- param.struc_p5.100
# +0.9 corr
Param.struc.1000[,,17] <- param.struc_p9.005
Param.struc.1000[,,18] <- param.struc_p9.010
Param.struc.1000[,,19] <- param.struc_p9.020
Param.struc.1000[,,20] <- param.struc_p9.100


#		SELECT
param.selec.1000 <- array(NA,c(1000,7,20))
# -0.9 corr
param.selec.1000[,,1] <- param.selec_m9.005
param.selec.1000[,,2] <- param.selec_m9.010
param.selec.1000[,,3] <- param.selec_m9.020
param.selec.1000[,,4] <- param.selec_m9.100
# -0.5 corr
param.selec.1000[,,5] <- param.selec_m5.005
param.selec.1000[,,6] <- param.selec_m5.010
param.selec.1000[,,7] <- param.selec_m5.020
param.selec.1000[,,8] <- param.selec_m5.100
# +0.0 corr
param.selec.1000[,,9] <- param.selec_p0.005
param.selec.1000[,,10] <- param.selec_p0.010
param.selec.1000[,,11] <- param.selec_p0.020
param.selec.1000[,,12] <- param.selec_p0.100
# +0.5 corr
param.selec.1000[,,13] <- param.selec_p5.005
param.selec.1000[,,14] <- param.selec_p5.010
param.selec.1000[,,15] <- param.selec_p5.020
param.selec.1000[,,16] <- param.selec_p5.100
# +0.9 corr
param.selec.1000[,,17] <- param.selec_p9.005
param.selec.1000[,,18] <- param.selec_p9.010
param.selec.1000[,,19] <- param.selec_p9.020
param.selec.1000[,,20] <- param.selec_p9.100



#		SELECTION & STRUCTURE
param.selStruc.1000 <- array(NA,c(1000,7,20))
# -0.9 corr
param.selStruc.1000[,,1] <- as.numeric(param.selStruc_m9.005)
param.selStruc.1000[,,2] <- as.numeric(param.selStruc_m9.010)
param.selStruc.1000[,,3] <- as.numeric(param.selStruc_m9.020)
param.selStruc.1000[,,4] <- as.numeric(param.selStruc_m9.100)
# -0.5 corr
param.selStruc.1000[,,5] <- as.numeric(param.selStruc_m5.005)
param.selStruc.1000[,,6] <- as.numeric(param.selStruc_m5.010)
param.selStruc.1000[,,7] <- as.numeric(param.selStruc_m5.020)
param.selStruc.1000[,,8] <- as.numeric(param.selStruc_m5.100)
# +0.0 corr
param.selStruc.1000[,,9] <- as.numeric(param.selStruc_p0.005)
param.selStruc.1000[,,10] <- as.numeric(param.selStruc_p0.010)
param.selStruc.1000[,,11] <- as.numeric(param.selStruc_p0.020)
param.selStruc.1000[,,12] <- as.numeric(param.selStruc_p0.100)
# +0.5 corr
param.selStruc.1000[,,13] <- as.numeric(param.selStruc_p5.005)
param.selStruc.1000[,,14] <- as.numeric(param.selStruc_p5.010)
param.selStruc.1000[,,15] <- as.numeric(param.selStruc_p5.020)
param.selStruc.1000[,,16] <- as.numeric(param.selStruc_p5.100)
# +0.9 corr
param.selStruc.1000[,,17] <- as.numeric(param.selStruc_p9.005)
param.selStruc.1000[,,18] <- as.numeric(param.selStruc_p9.010)
param.selStruc.1000[,,19] <- as.numeric(param.selStruc_p9.020)
param.selStruc.1000[,,20] <- as.numeric(param.selStruc_p9.100)



#	Correlation 

correlation.1000 <- cbind(correlation.005,correlation.010,correlation.020,correlation.100)


save(Param.struc.1000, param.selec.1000,param.selStruc.1000,correlation.1000,file="Data_Sample_1000.RData")



#### 


####


#### 


####

#### 

v


# P` L ` O ` T ` T ` I ` N ` G

# we are using this for plot: ###### XXXX_-XX_010

rm(list=ls())
load("/Users/lucas/Documents/Academia/Columbia University/Methods Paper/Selection/Results and Plots/Data_Sample_1000.RData")

objects()

par(family="CMU Serif", mfrow=c(1,2))
#par(family="CMU Serif")
plot(mean(correlation.1000[,6]), mean(as.numeric(Param.struc.1000[,6,2])-1), xlim=c(-1,1), ylim=c(-0.05,0.35), pch=16, col=rgb(255,165,0,250,maxColorValue=255), cex=2.6, ylab="Bias", xlab="Correlation of Errors", main=expression(paste("Estimate of ", beta[21])))
points(mean(correlation.1000[,6]),mean(as.numeric(param.selec.1000[,6,2])-1),pch=16, col=rgb(0,139,139,200,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,6]),mean(as.numeric(param.selStruc.1000[,6,2])-1),pch=16, col=rgb(160,32,240,200,maxColorValue=255), cex=2.6)
abline(h=0, lty=2)

points(mean(correlation.1000[,7]), mean(as.numeric(Param.struc.1000[,6,6])-1), pch=16, col=rgb(255,165,0,250,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,7]),mean(as.numeric(param.selec.1000[,6,6])-1),pch=16, col=rgb(0,139,139,200,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,7]),mean(as.numeric(param.selStruc.1000[,6,6])-1),pch=16, col=rgb(160,32,240,200,maxColorValue=255), cex=2.6)


points(mean(correlation.1000[,8]), mean(as.numeric(Param.struc.1000[,6,10])-1), pch=16, col=rgb(255,165,0,250,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,8]),mean(as.numeric(param.selec.1000[,6,10])-1),pch=16, col=rgb(0,139,139,200,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,8]),mean(as.numeric(param.selStruc.1000[,6,10])-1),pch=16, col=rgb(160,32,240,200,maxColorValue=255), cex=2.6)


points(mean(correlation.1000[,9]), mean(as.numeric(Param.struc.1000[,6,14])-1), pch=16, col=rgb(255,165,0,250,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,9]),mean(as.numeric(param.selec.1000[,6,14])-1),pch=16, col=rgb(0,139,139,200,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,9]),mean(as.numeric(param.selStruc.1000[,6,14])-1),pch=16, col=rgb(160,32,240,200,maxColorValue=255), cex=2.6)


points(mean(correlation.1000[,10]), mean(as.numeric(Param.struc.1000[,6,18])-1), pch=16, col=rgb(255,165,0,250,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,10]),mean(as.numeric(param.selec.1000[,6,18])-1),pch=16, col=rgb(0,139,139,200,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,10]),mean(as.numeric(param.selStruc.1000[,6,18])-1),pch=16, col=rgb(160,32,240,200,maxColorValue=255), cex=2.6)

#	#	#	#	#	#	#	#	#	#	#	#	#	#	#
# Lines for the plot
# yellow
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(Param.struc.1000[,6,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(Param.struc.1000[,6,6])-1), lty=2, col=rgb(255,165,0,250,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(Param.struc.1000[,6,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(Param.struc.1000[,6,10])-1), lty=2, col=rgb(255,165,0,250,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(Param.struc.1000[,6,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(Param.struc.1000[,6,14])-1), lty=2, col=rgb(255,165,0,250,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(Param.struc.1000[,6,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(Param.struc.1000[,6,18])-1), lty=2, col=rgb(255,165,0,250,maxColorValue=255), lwd=2)

# green
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(param.selec.1000[,6,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(param.selec.1000[,6,6])-1), lty=2, col=rgb(0,139,139,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(param.selec.1000[,6,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(param.selec.1000[,6,10])-1), lty=2, col=rgb(0,139,139,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(param.selec.1000[,6,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(param.selec.1000[,6,14])-1), lty=2, col=rgb(0,139,139,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(param.selec.1000[,6,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(param.selec.1000[,6,18])-1), lty=2, col=rgb(0,139,139,200,maxColorValue=255), lwd=2)

# purple
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(param.selStruc.1000[,6,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(param.selStruc.1000[,6,6])-1), lty=2, col=rgb(160,32,240,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(param.selStruc.1000[,6,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(param.selStruc.1000[,6,10])-1), lty=2, col=rgb(160,32,240,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(param.selStruc.1000[,6,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(param.selStruc.1000[,6,14])-1), lty=2, col=rgb(160,32,240,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(param.selStruc.1000[,6,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(param.selStruc.1000[,6,18])-1), lty=2, col=rgb(160,32,240,200,maxColorValue=255), lwd=2)


legend(.4,.3, legend=c("Stc", "HSl", "StSl"), col=c(rgb(255,165,0,250,maxColorValue=255),rgb(0,139,139,200,maxColorValue=255),rgb(160,32,240,200,maxColorValue=255)), cex=1.2, pch=16, bty="n")


### 2nd plot

plot(mean(correlation.1000[,6]), mean(as.numeric(Param.struc.1000[,5,2])-1), xlim=c(-1,1), ylim=c(-0.05,0.35), pch=16, col=rgb(255,165,0,250,maxColorValue=255), cex=2.6, ylab="Bias", xlab="Correlation of Errors", main=expression(paste("Estimate of ", beta[22])))
points(mean(correlation.1000[,6]),mean(as.numeric(param.selec.1000[,5,2])-1),pch=16, col=rgb(0,139,139,200,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,6]),mean(as.numeric(param.selStruc.1000[,5,2])-1),pch=16, col=rgb(160,32,240,200,maxColorValue=255), cex=2.6)
abline(h=0, lty=2)

points(mean(correlation.1000[,7]), mean(as.numeric(Param.struc.1000[,5,6])-1), pch=16, col=rgb(255,165,0,250,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,7]),mean(as.numeric(param.selec.1000[,5,6])-1),pch=16, col=rgb(0,139,139,200,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,7]),mean(as.numeric(param.selStruc.1000[,5,6])-1),pch=16, col=rgb(160,32,240,200,maxColorValue=255), cex=2.6)


points(mean(correlation.1000[,8]), mean(as.numeric(Param.struc.1000[,5,10])-1), pch=16, col=rgb(255,165,0,250,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,8]),mean(as.numeric(param.selec.1000[,5,10])-1),pch=16, col=rgb(0,139,139,200,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,8]),mean(as.numeric(param.selStruc.1000[,5,10])-1),pch=16, col=rgb(160,32,240,200,maxColorValue=255), cex=2.6)


points(mean(correlation.1000[,9]), mean(as.numeric(Param.struc.1000[,5,14])-1), pch=16, col=rgb(255,165,0,250,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,9]),mean(as.numeric(param.selec.1000[,5,14])-1),pch=16, col=rgb(0,139,139,200,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,9]),mean(as.numeric(param.selStruc.1000[,5,14])-1),pch=16, col=rgb(160,32,240,200,maxColorValue=255), cex=2.6)


points(mean(correlation.1000[,10]), mean(as.numeric(Param.struc.1000[,5,18])-1), pch=16, col=rgb(255,165,0,250,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,10]),mean(as.numeric(param.selec.1000[,5,18])-1),pch=16, col=rgb(0,139,139,200,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,10]),mean(as.numeric(param.selStruc.1000[,5,18])-1),pch=16, col=rgb(160,32,240,200,maxColorValue=255), cex=2.6)

#	#	#	#	#	#	#	#	#	#	#	#	#	#	#
# Lines for the plot
# yellow
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(Param.struc.1000[,5,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(Param.struc.1000[,5,6])-1), lty=2, col=rgb(255,165,0,250,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(Param.struc.1000[,5,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(Param.struc.1000[,5,10])-1), lty=2, col=rgb(255,165,0,250,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(Param.struc.1000[,5,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(Param.struc.1000[,5,14])-1), lty=2, col=rgb(255,165,0,250,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(Param.struc.1000[,5,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(Param.struc.1000[,5,18])-1), lty=2, col=rgb(255,165,0,250,maxColorValue=255), lwd=2)

# green
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(param.selec.1000[,5,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(param.selec.1000[,5,6])-1), lty=2, col=rgb(0,139,139,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(param.selec.1000[,5,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(param.selec.1000[,5,10])-1), lty=2, col=rgb(0,139,139,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(param.selec.1000[,5,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(param.selec.1000[,5,14])-1), lty=2, col=rgb(0,139,139,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(param.selec.1000[,5,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(param.selec.1000[,5,18])-1), lty=2, col=rgb(0,139,139,200,maxColorValue=255), lwd=2)

# purple
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(param.selStruc.1000[,5,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(param.selStruc.1000[,5,6])-1), lty=2, col=rgb(160,32,240,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(param.selStruc.1000[,5,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(param.selStruc.1000[,5,10])-1), lty=2, col=rgb(160,32,240,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(param.selStruc.1000[,5,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(param.selStruc.1000[,5,14])-1), lty=2, col=rgb(160,32,240,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(param.selStruc.1000[,5,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(param.selStruc.1000[,5,18])-1), lty=2, col=rgb(160,32,240,200,maxColorValue=255), lwd=2)




legend(.4,.3, legend=c("Stc", "HSl", "StSl"), col=c(rgb(255,165,0,250,maxColorValue=255),rgb(0,139,139,200,maxColorValue=255),rgb(160,32,240,200,maxColorValue=255)), cex=1.2, pch=16, bty="n")
####

#### 


####

#### 


####

#par(family="CMU Serif", mfrow=c(1,2))
par(family="CMU Serif")
plot(mean(correlation.1000[,6]), mean(as.numeric(Param.struc.1000[,3,2])-1), xlim=c(-1,1), ylim=c(-0.7,0.4), pch=16, col=rgb(255,165,0,250,maxColorValue=255), cex=2.6, ylab="Bias", xlab="Correlation of Errors", main=expression(paste("Estimate of ", beta[11])))
points(mean(correlation.1000[,6]),mean(as.numeric(param.selec.1000[,3,2])-1),pch=16, col=rgb(0,139,139,200,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,6]),mean(as.numeric(param.selStruc.1000[,3,2])-1),pch=16, col=rgb(160,32,240,200,maxColorValue=255), cex=2.6)
abline(h=0, lty=2)

points(mean(correlation.1000[,7]), mean(as.numeric(Param.struc.1000[,3,6])-1), pch=16, col=rgb(255,165,0,250,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,7]),mean(as.numeric(param.selec.1000[,3,6])-1),pch=16, col=rgb(0,139,139,200,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,7]),mean(as.numeric(param.selStruc.1000[,3,6])-1),pch=16, col=rgb(160,32,240,200,maxColorValue=255), cex=2.6)


points(mean(correlation.1000[,8]), mean(as.numeric(Param.struc.1000[,3,10])-1), pch=16, col=rgb(255,165,0,250,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,8]),mean(as.numeric(param.selec.1000[,3,10])-1),pch=16, col=rgb(0,139,139,200,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,8]),mean(as.numeric(param.selStruc.1000[,3,10])-1),pch=16, col=rgb(160,32,240,200,maxColorValue=255), cex=2.6)


points(mean(correlation.1000[,9]), mean(as.numeric(Param.struc.1000[,3,14])-1), pch=16, col=rgb(255,165,0,250,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,9]),mean(as.numeric(param.selec.1000[,3,14])-1),pch=16, col=rgb(0,139,139,200,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,9]),mean(as.numeric(param.selStruc.1000[,3,14])-1),pch=16, col=rgb(160,32,240,200,maxColorValue=255), cex=2.6)


points(mean(correlation.1000[,10]), mean(as.numeric(Param.struc.1000[,3,18])-1), pch=16, col=rgb(255,165,0,250,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,10]),mean(as.numeric(param.selec.1000[,3,18])-1),pch=16, col=rgb(0,139,139,200,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,10]),mean(as.numeric(param.selStruc.1000[,3,18])-1),pch=16, col=rgb(160,32,240,200,maxColorValue=255), cex=2.6)

#	#	#	#	#	#	#	#	#	#	#	#	#	#	#
# Lines for the plot
# yellow
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(Param.struc.1000[,3,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(Param.struc.1000[,3,6])-1), lty=2, col=rgb(255,165,0,250,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(Param.struc.1000[,3,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(Param.struc.1000[,3,10])-1), lty=2, col=rgb(255,165,0,250,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(Param.struc.1000[,3,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(Param.struc.1000[,3,14])-1), lty=2, col=rgb(255,165,0,250,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(Param.struc.1000[,3,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(Param.struc.1000[,3,18])-1), lty=2, col=rgb(255,165,0,250,maxColorValue=255), lwd=2)

# green
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(param.selec.1000[,3,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(param.selec.1000[,3,6])-1), lty=2, col=rgb(0,139,139,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(param.selec.1000[,3,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(param.selec.1000[,3,10])-1), lty=2, col=rgb(0,139,139,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(param.selec.1000[,3,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(param.selec.1000[,3,14])-1), lty=2, col=rgb(0,139,139,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(param.selec.1000[,3,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(param.selec.1000[,3,18])-1), lty=2, col=rgb(0,139,139,200,maxColorValue=255), lwd=2)

# purple
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(param.selStruc.1000[,3,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(param.selStruc.1000[,3,6])-1), lty=2, col=rgb(160,32,240,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(param.selStruc.1000[,3,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(param.selStruc.1000[,3,10])-1), lty=2, col=rgb(160,32,240,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(param.selStruc.1000[,3,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(param.selStruc.1000[,3,14])-1), lty=2, col=rgb(160,32,240,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(param.selStruc.1000[,3,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(param.selStruc.1000[,3,18])-1), lty=2, col=rgb(160,32,240,200,maxColorValue=255), lwd=2)



legend(.4,.4, legend=c("Stc", "HSl", "StSl"), col=c(rgb(255,165,0,250,maxColorValue=255),rgb(0,139,139,200,maxColorValue=255),rgb(160,32,240,200,maxColorValue=255)), cex=1.2, pch=16, bty="n")

### 2nd plot

plot(mean(correlation.1000[,6]), mean(as.numeric(Param.struc.1000[,1,2])-1), xlim=c(-1,1), ylim=c(-2.1,1.4), pch=16, col=rgb(255,165,0,250,maxColorValue=255), cex=2.6, ylab="Bias", xlab="Correlation of Errors",  main=expression(paste("Estimate of ", beta[12])))
points(mean(correlation.1000[,6]),mean(as.numeric(param.selec.1000[,2,2])-1),pch=16, col=rgb(0,139,139,200,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,6]),mean(as.numeric(param.selStruc.1000[,1,2])-1),pch=16, col=rgb(160,32,240,200,maxColorValue=255), cex=2.6)
abline(h=0, lty=2)

points(mean(correlation.1000[,7]), mean(as.numeric(Param.struc.1000[,1,6])-1), pch=16, col=rgb(255,165,0,250,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,7]),mean(as.numeric(param.selec.1000[,2,6])-1),pch=16, col=rgb(0,139,139,200,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,7]),mean(as.numeric(param.selStruc.1000[,1,6])-1),pch=16, col=rgb(160,32,240,200,maxColorValue=255), cex=2.6)


points(mean(correlation.1000[,8]), mean(as.numeric(Param.struc.1000[,1,10])-1), pch=16, col=rgb(255,165,0,250,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,8]),mean(as.numeric(param.selec.1000[,2,10])-1),pch=16, col=rgb(0,139,139,200,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,8]),mean(as.numeric(param.selStruc.1000[,1,10])-1),pch=16, col=rgb(160,32,240,200,maxColorValue=255), cex=2.6)


points(mean(correlation.1000[,9]), mean(as.numeric(Param.struc.1000[,1,14])-1), pch=16, col=rgb(255,165,0,250,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,9]),mean(as.numeric(param.selec.1000[,2,14])-1),pch=16, col=rgb(0,139,139,200,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,9]),mean(as.numeric(param.selStruc.1000[,1,14])-1),pch=16, col=rgb(160,32,240,200,maxColorValue=255), cex=2.6)


points(mean(correlation.1000[,10]), mean(as.numeric(Param.struc.1000[,1,18])-1), pch=16, col=rgb(255,165,0,250,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,10]),mean(as.numeric(param.selec.1000[,2,18])-1),pch=16, col=rgb(0,139,139,200,maxColorValue=255), cex=2.6)
points(mean(correlation.1000[,10]),mean(as.numeric(param.selStruc.1000[,1,18])-1),pch=16, col=rgb(160,32,240,200,maxColorValue=255), cex=2.6)

#	#	#	#	#	#	#	#	#	#	#	#	#	#	#
# Lines for the plot
# yellow
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(Param.struc.1000[,1,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(Param.struc.1000[,1,6])-1), lty=2, col=rgb(255,165,0,250,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(Param.struc.1000[,1,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(Param.struc.1000[,1,10])-1), lty=2, col=rgb(255,165,0,250,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(Param.struc.1000[,1,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(Param.struc.1000[,1,14])-1), lty=2, col=rgb(255,165,0,250,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(Param.struc.1000[,1,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(Param.struc.1000[,1,18])-1), lty=2, col=rgb(255,165,0,250,maxColorValue=255), lwd=2)

# green
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(param.selec.1000[,2,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(param.selec.1000[,2,6])-1), lty=2, col=rgb(0,139,139,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(param.selec.1000[,2,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(param.selec.1000[,2,10])-1), lty=2, col=rgb(0,139,139,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(param.selec.1000[,2,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(param.selec.1000[,2,14])-1), lty=2, col=rgb(0,139,139,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(param.selec.1000[,2,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(param.selec.1000[,2,18])-1), lty=2, col=rgb(0,139,139,200,maxColorValue=255), lwd=2)

# purple
segments(x0=mean(correlation.1000[,6]), y0=mean(as.numeric(param.selStruc.1000[,1,2])-1), x=mean(correlation.1000[,7]), y=mean(as.numeric(param.selStruc.1000[,1,6])-1), lty=2, col=rgb(160,32,240,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,7]), y0=mean(as.numeric(param.selStruc.1000[,1,6])-1), x=mean(correlation.1000[,8]), y=mean(as.numeric(param.selStruc.1000[,1,10])-1), lty=2, col=rgb(160,32,240,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,8]), y0=mean(as.numeric(param.selStruc.1000[,1,10])-1), x=mean(correlation.1000[,9]), y=mean(as.numeric(param.selStruc.1000[,1,14])-1), lty=2, col=rgb(160,32,240,200,maxColorValue=255), lwd=2)

segments(x0=mean(correlation.1000[,9]), y0=mean(as.numeric(param.selStruc.1000[,1,14])-1), x=mean(correlation.1000[,10]), y=mean(as.numeric(param.selStruc.1000[,1,18])-1), lty=2, col=rgb(160,32,240,200,maxColorValue=255), lwd=2)



legend(.4,1.4, legend=c("Stc", "HSl", "StSl"), col=c(rgb(255,165,0,250,maxColorValue=255),rgb(0,139,139,200,maxColorValue=255),rgb(160,32,240,200,maxColorValue=255)), cex=1.2, pch=16, bty="n")


#### Correlation:


# sample 200
load("/Users/lucas/Documents/Academia/Columbia University/Methods Paper/Selection/Results and Plots/Data_Sample_0200.RData")
coRR <- c(correlation.0200[,6],correlation.0200[,7],correlation.0200[,8],correlation.0200[,9],correlation.0200[,10])
hist(coRR, breaks=200, col=rgb(130,0,255,120,maxColorValue=255), border=FALSE, ylim=c(0,550))#, add=TRUE)



# sample 500
load("/Users/lucas/Documents/Academia/Columbia University/Methods Paper/Selection/Results and Plots/Data_Sample_0500.RData")
coRR <- c(correlation.0500[,6],correlation.0500[,7],correlation.0500[,8],correlation.0500[,9],correlation.0500[,10])
hist(coRR, breaks=200, col=rgb(130,0,255,120,maxColorValue=255), border=FALSE, add=TRUE)



load("/Users/lucas/Documents/Academia/Columbia University/Methods Paper/Selection/Results and Plots/Data_Sample_1000.RData")
coRR <- c(correlation.1000[,6],correlation.1000[,7],correlation.1000[,8],correlation.1000[,9],correlation.1000[,10])
hist(coRR, breaks=200, col=rgb(130,0,255,120,maxColorValue=255), border=FALSE, add=TRUE)



rho.m9 <- as.numeric(param.selStruc.1000[,7,2])
rho.m5 <- as.numeric(param.selStruc.1000[,7,6])
rho.p0 <- as.numeric(param.selStruc.1000[,7,10])
rho.p5 <- as.numeric(param.selStruc.1000[,7,14])
rho.p9 <- as.numeric(param.selStruc.1000[,7,18])


hist(rho.m9, breaks=100, xlim=c(-1,1), ylim=c(0,200),border=FALSE, col=rgb(0,0,255,100,maxColorValue=255))
hist(rho.m5, breaks=100, add=TRUE, border=FALSE, col=rgb(0,0,255,100,maxColorValue=255))
hist(rho.p0, breaks=100, add=TRUE, border=FALSE, col=rgb(0,0,255,100,maxColorValue=255))
hist(rho.p5, breaks=100, add=TRUE, border=FALSE, col=rgb(0,0,255,100,maxColorValue=255))
hist(rho.p9, breaks=100, add=TRUE, border=FALSE, col=rgb(0,0,255,100,maxColorValue=255))
rHO <- c(rho.m9,rho.m5,rho.p0,rho.p5,rho.p9)
hist(coRR, breaks=200, col=rgb(130,0,255,120,maxColorValue=255), border=FALSE, add=TRUE)

par(family="CMU Serif")
hist(coRR, breaks=300, col=rgb(0,0,255,180,maxColorValue=255), border=FALSE, xlim=c(-1,1), main=expression(paste("Estimate of ", rho," based on AESe")), xlab="sdsdsdsd")
hist(rHO, breaks=400, add=TRUE)
